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Abstract 

We consider the 2D J\ — J2 classical XY model on a square lattice. In the 
frustrated phase corresponding to J2 > J\/2, an Ising like order parameter 
emerges by an "order due to disorder" effect. This leads to a discrete Z2 
symmetry plus the U(l) global one. Using a powerful algorithm we show 
that the system undergoes two transitions at different but still very close 
temperatures, one of Kosterlitz-Thouless (KT) type and another one which 
does not belong to the expected Ising universality class. A new analysis 
of the KT transition has been developed in order to avoid the use of the 
non-universal helicity jump and to allow the computation of the exponents 
without a precise determination of the critical temperature. Moreover, our 
huge number of data enables us to exhibit the existence of large finite size 
effects explaining the dispersed results found in the literature concerning the 
more studied frustrated 2D, XY models. 
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I. INTRODUCTION 



The ground state of a large class of two-dimensional classical frustrated XY models have 
the particularity to exhibit both continuous and discrete Z 2 degeneracy simultaneously in the 
ground state. It results the appearance of a new Ising-like order parameter, in addition to the 
continuous U(l) symmetry. The most famous example exhibiting such behaviors is certainly 
the fully frustrated XY model (FFXY) which was originally introduced by Villain as a 
frustrated XY model without disorder PJ . In this model, the Z2 symmetry is associated with 
the two types of chirality ordering. This model is also of great interest because it describes 
a superconducting array of Josephson junctions under an external transverse magnetic field 
such that the flux per plaquette is half the quantum flux f2|. For fifteen years, extensive 



(essentially numerical) works have been carried on the FFXY ||- [ 10 1 and also to some 
related models like the triangular antiferromagnetic XY model [[11]], the helical XY model 



12], the Coulomb gas system of half integer charges fl3|,|14|. The interplay between the two 



transitions can lead a priori to two transitions, namely a Kosterlitz-Thouless one and an 
Ising one. Nevertheless, the entanglement between the two order parameters considerably 
complicates the analysis. The nature of the phase diagram is still rather unconclusive and 
controversial. Three different scenari have been advocated: either the two transitions occur 
at the same point and eventually merge to give a new universality class or the 

two transitions occur at different points and are of Ising and Kosterlitz-Thouless types plus 
some strong finite size effects |]13 , |15| ; or finally the two transitions are effectively separated 
but the transition associated to the chiral order parameter is not of Ising type [^,pj|l4| . The 
most recent numerical studies are in favor of the latter scenario. Nonetheless, without a 
strong analytical support, the problem is still completely open. 

The purpose of the present article is other. We want to clarify the critical behavior of 
a less studied frustrated XY model: the 2D, J\ — J2 XY model on a square lattice which, 
as will be shown, is in the same universality class as the models quoted above (or more 
precisely has the same problematic). The Hamiltonian reads 
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h = -j 1 J2 + ^2 E s ^ (i) 

<i,j> «h,l» 

= -J, £ cos(^ - 0,) + J 2 £ cos (^ - °i) » ( 2 ) 

<ij> «k,l» 

where S, are two component classical vectors of unit length, with Ji, J 2 > 0, < > and 
<< >> indicate respectively the sum over nearest neighbors and next to nearest neighbors. 
When 2Ji > J 2 the ground state is ferromagnetic. It leads to a Kosterlitz-Thouless (KT) 
transition at the temperature T^t ~ n ( Jl ~ 2j2 ^ [|j~6[| . This temperature is obtained from 
the Villain treatment of the Hamiltonian (JJ) (see ref. |l7j for details). However, when 
2Ji < J 2 , the ground state consists of two independent \[2 x \/2 sublattices with AF order. 
The ground state energy E = — 2iVJ 2 does not depend on 0, an angle parameterizing 
the relative orientations between both sublattices. This non-trivial degeneracy is lifted by 
thermal fluctuations and a collinear ordering (corresponding to = or it) is selected 
|18| . The two possible ground states are depicted in Fig. [1|. The angle thus plays a 
role analogous to the chiral order parameter. This selection mechanism is one of the most 



famous "order due to disorder" effect |18[ in a sense that fluctuations brings kind of order 
by lifting this extra continuous symmetry. The resulting symmetry is therefore U(l) x Z 2 . 
Monte-Carlo simulations predict only a low temperature phase with a nematic ordering and 
a disordered high temperature phase [flS|, 19(1 - The critical behavior has, as far as we know, 



only been partially explored in ref. [|l9j. Unfortunately the results are very approximate 
and no definite conclusion on the presence of one or two transitions and either on their 
universality classes has been given. In this work, we have carried on extensive numerical 
Monte Carlo simulations on the J\ — J 2 XY model using algorithms which allow us to obtain 
very accurate and robust results. 

We now give the outline of the paper. Section 2 contains a brief summary of analytical 
results showing the relations between the J\ — J 2 XY model and the Ising-XY model, which 
is a generic model used to describe the universality class of frustrated XY models with 
symmetry Z 2 x U(l) like the FFXY. In section 3, we present our numerical results and the 
analysis of some critical exponents. Finally, section 4 is devoted to the discussion of the 



results and to a brief conclusion. The estimation of statistical and systematic errors has 
been relegated in the appendices. 



II. THE Ji - J 2 XY AND ISING-XY MODEL 

In this section, we sum up the main analytical results concerning the J\ — J 2 XY model. 
We essentially focus on the more interesting frustrated phase corresponding to 2J 2 > J%, 
where the ground state consists of two independent sublattices. Thermal fluctuations select 
a collinear ordering |18 |, and we have two kinds of domains represented in Fig. |l| The first 



step, following Chandra et al. ||20|| , is to perform a gradient expansion of the classical energy 
(H|). The problem is now translated in a new one on a (2 x 2) square lattice, but now with 
two spins 1 and 2 per vertices pointing in the same directions. The new classical action A 
reads 

^ = -^E[W) 2 + (w 2 ) 2 
+ 2Acos0 (v x e 1 v x e 2 -v y e 1 v y e 2 )} , (3) 



where we have defined A = ^j- < 1 and introduced the lattice derivatives V x , V y |L6 



The signature of the U(l) degeneracy lies now in the strong anisotropy between x and y 
directions. The cos0 labels the different possible classical ground states at T = 0. Notice 
that, if we do the Gaussian integration over the angular variables, we recover the results of 
Henley |18|| , namely 

A ~ const- 0.32 (A^V , (4) 
\ 2 J 2 ) 

proving that a collinear ordering is selected when minimizing the free energy(cos(0) = ±1). 

Let us now include the effects of the vortices. The most natural way to include them 
would be to apply the Villain transformation to all quadratic terms in the action (||). Such 
a treatment is quite tedious and unappropriated because the vortices variables built with 
the anisotropic term (V x 6 l 1 V :r 6 l 2 — V y 6' 1 V y 6 l 2 ) are not well defined due to the strong fluc- 
tuations between the two sublattices which tend to decouple in the infrared limit ( see 



16] for more details). A simplest way to take into account the coupling between the two 
sublattices is to replace the anisotropic term in ([$]) by the local spin waves effect i.e. by 
— 0.32A 2 cos 2 (#i(r) — 6 2 (r)). This term is just the local version of (0). Such a treatment 



has already been used by Garel et al. for helimagnets [JT2J] and also by Chandra et al. for 
the J% — J 2 Heisenberg model [ZIJ. By applying the Villain transformation to the first two 
terms in @ and using 



exp [hcosp(9i(r) — O2 if))} = E ex P [ipS(r)(0i(f) — 62(f)) 

S(r) 

+ \ogy s S 2 (r) , 



(5) 



with p = 2 and y s = h/2 = 0.08A 2 , we obtain 



Z= E EE /W^expf-^EE [(V^(r)-2vr<(r)) 2 

{n5*(r)}K(r)}S(r) J V r i=1 > 2 

+ E [^W^iW - 2 (r)) + logy s 5 2 (f) 



(6) 



The nf (z = 1, 2) are integer link variables living on the two diagonal sublattices. By 
integrating on angular variables, we easily find 



A//= E 



n(3J 2 M l {r) log ^— ^M^r') + 7i(3J 2 M 2 {r) log Jl_^M 2 (r') 



2 1 'i 

-ip{M x {r) +M 2 (r))Q\r-r'\S(r') + ^S(r) log lT ~ rl S(r') 

27rp a 

+ E [log Vi (M l f(r) + \ogy2 (M 2 ) 2 (r) + \ogy s S 2 {r) 



(?) 



where we have introduced the vortex variables Mi = e^V^n" corresponding to vortices 
on the two sublattices. The fugacities are as usual considered as genuine variables defined 



initially by = exp( 



. The interaction is defined by 0|r — r'\ = arctan(^3-), and 



verifies d y log |r — r'\ = — d x Q\r — r'\. This action corresponds to two coupled XY models. 
Under real space renormalization group transformations, the coupling term is strongly rele- 



vant and locks the phase difference in 6\{r) = #2(7") + kir with k — 0, 1 [21|. It leads in the 
strong coupling limit to an effective model whose Hamiltonian has the following form 
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Hi-xy = - £ [A(l + a^j) cos(^ - 6,) + Cwj] . (8) 

<*J> 

The value of A and C depend of the initial values of h = f(J±, J2) and (3J2- This model 
refereed as the Ising-XY model in the literature has been largely debated. This Ising-XY 
model is believed to describe the critical behavior of all frustrated XY models quoted in the 
introduction. The phase diagram has been obtained numerically by Granato et al. and 
has been reproduced for convenience in Fig. ^. Three different phases can be distinguished: 
the upper right corner phase correspond to the low temperature ordered phase, the low 
left corner phase is the high temperature disordered one, and the intermediate one is Ising 
ordered but XY disordered (namely with free vortices). Above the point P, the Ising and XY 
transitions are well separated and mix under P. The question concerning the transition(s) 
under P is still under debate. A recent work of S. Lee et al. seems to indicate that 
the two transitions never merge completely but get closer [15||. Nevertheless the Ising-like 
magnetization exponent has been found different from 1 and continuously varying along the 
line (PT) |[22|| . We have shown that the J\ — J2 XY model should be also described by the 
Ising-XY model and should therefore correspond to a curve crossing the line under P (so 
with only one or two very close transitions). Since we are not able to provide analytical 
relations between (h, J2) and (A,C), the form of this curve and its intersection with the 
segment (PT) is unaccessible. Moreover, when varying J2/J1, we shall obtain a different 



intersection point as it was firstly noted in |jl6 |. Nevertheless, it opens the possibility that 
the critical exponent v should vary with the ratio J2/J1 as in the analysis of Granato al. 



22] or of Lee al. |15|. Similar considerations have been done in the study of a generalized 



frustrated XY model where an extra-parameter has been added ||10|| . No clear conclusion 
concerning the nature of the phase transitions can be therefore drawn at this level. The 
purpose of the next section is therefore to answer these questions by help of extensive Monte 
Carlo simulations. Moreover, it can also be regarded as an indirect way of studying the 
Ising-XY model and other related models. 
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III. MONTE CARLO ANALYSIS 



A. Observables 

As explained above we can define two order parameters corresponding to the two sym- 
metries U(l) and Z2. The first one is the total magnetization M defined by the sum of all 
spins, the second is the chirality k defined by the sum of all chiralities Ki defined on each 
cell by: 

Kt = ±(Si-S k ){Sj-Si) (9) 

where (i, j, k, I) are the four corners of one cell with diagonal (i,k) and (j,l). The two ground 
states depicted in Fig. [l] have Ki — ±1. 

We have studied our system in the finite size scaling region where the correlation length 
is much bigger than the lattice size. The quantities needed for our analysis are defined 
below. For each temperature we calculate: 

^ (<K V K>2) 

iV<^> 

k B T ( ' 

V," = i^i- < E > (13) 

< K > 

V« = <^> - <E > (15) 

U» = l- <M4> (16) 
3 < M 2 > 2 1 ; 

U K = 1 - <K > (17) 

where E is the energy, x is the magnetic susceptibility per site, V\$ are cumulants used to 
obtain the critical exponent u, U are the fourth order cumulants, < ... > means the thermal 
average. 



B. Algorithm 



We use in this work the standard Metropolis algorithm. At each site a new random 
orientation for the spin is chosen. The interaction energy between this spin and its neighbors 
is then calculated. If lower than the energy of the old state, the new state is accepted, 
otherwise, it is accepted only with a probability p according to the standard Metropolis 
algorithm. 

However the critical slowing down is important and we improve the speed of the simula- 
tion using the local over- relaxation algorithm (OR) |23|| . This algorithm consists in choosing 
the new orientation of the spin such that the energy remains unchanged. For XY spins the 
only possibility is to take the symmetric of the old spin to the local field (the sum of the 
neighbors). This algorithm is obviously non ergodic, i.e. all states can not be reached. It 
must thus be used in combination with the standard Metropolis algorithm (MET). There- 
fore, at each step (regarded as one unit) we use one MET step and a certain number of steps 
of over-relaxation (NOR) algorithm. The larger NOR, the smaller the autocorrelation time 
(the number of steps between two independent configurations), but then the larger the time 
needed for each step is. We have thus to choose the best NOR to minimize the real autocor- 
relation time. This depends on the time needed for each algorithm. In our implementation 
the Metropolis algorithm is six times larger than the over relaxation algorithm. 

In order to calculate the autocorrelation time we follow the procedure explained in ap- 
pendix [A]. In table | we present the results of the autocorrelation time r K for different NOR 
at the critical temperature T C K for a lattice size L = 30 and J2/ J\ = 0.7. The second column 
gives t k while the third column gives the real autocorrelation time t k (1 + NOR/6), i.e. 
the quantity to be minimized. The value NOR ~ L/5 = 6 seems to fit best. We have 
checked that this ratio does not change significantly for sizes L = 20 and L = 40, which is 
in agreement with the argument of Adler |24j] where the best NOR should be proportional 
to the correlation length, i.e. to the lattice size in the finite size scaling region. 

In Fig |3| we have shown in a log-log scale the real chirality autocorrelation time function 
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of the lattice size for NOR = and NOR = L/5. For larger lattice sizes the gain is more 
than a factor 30 using the over relaxation algorithm. The critical exponent z defined by 
t ~ L z is 2.29(4) without the use of the OR algorithm and is in agreement with the results 



on the Villain lattice 2.31 but in disagreement with the dynamic approach of Luo et al. 
H who obtained 2.17(4). We note for this last case that an error on z leads to errors on the 
other exponents. 

In the following the simulations have been done using NOR = L/5 for each Metropolis 
algorithm. For each simulation, we use a number t av measurements, made after an updating 
time t up is carried out for equilibration. For each size, between three and six different 
initial configurations (ordered or random) have been tested to be sure that our system is 
not trapped into a metastable state. In table H we present some details of our simulations. 
We want to stress that the number of Monte Carlo steps used in this work is one order of 
magnitude larger than previous studies and, combined with a better algorithm, produces a 
better estimate of the thermodynamic quantities. 

Our errors are calculated with the help of the Jackknife procedure [^BJ. When compiling 
the different results from previous studies we have noticed that the errors reported are quite 
often strongly underestimated. Therefore we have presented in the appendix ^ our method 
to evaluate the errors coming from the simulation and in particular a simplified method of 
the Jackknife procedure. 

We use in this work the histogram MC technique developed by Ferrenberg and Swendsen 



27,28]. From a simulation done at T , this technique allows to obtain thermodynamic 



quantities at T close to T . 



C. Phase diagram 

We have performed many simulations in varying the value of J 2/ J\ to obtain the critical 
temperature T c which is represented in Fig [|. The transition for J2/J1 < 0.5 is a standard 
Kosterlitz-Thouless transition in agreement with theoretical predictions. For J2/J1 > 0.5 it 
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is difficult to discriminate between the hypothesis of one or two transitions separating the low 
temperature nematic phase from the high temperature disordered phase. We have therefore 
decided to focus on the particular value J2 = 0.7, Ji = 1. (black circle) in the remainder 
of this work. It is worth noticing that Fernandez et al. |B| have done their calculation for 
J 2 = Ji = 1. 

As it was first emphasized in section 2, it is possible that the exponents could vary with 
the ratio J2/J1 fL6|. This should be coherent with the picture proposed by Granato et al for 



the Ising-Xy model |22j and by Lee et al. [15| . Nevertheless, the first step is to perform very 



highly accurate Monte Carlo simulations at some fixed value of J2/J1 to show the existence 
to test the existence of two close transitions, and check that the chiral magnetic exponent v 
is clearly different from 1. 

D. Z2 symmetry 

We concentrate first on the breakdown of the Z2 symmetry with the order parameter k 
defined in (|). 

To find the critical temperature T c we record the variation of U K with T for various 
system sizes in Fig |] and then locate T C K at the intersection of these curves |29| since the 
ratio of U K for two different lattice sizes L and V = bL should be 1 at T C K , namely 



TTK 
U bL 



= 1. (18) 

T=T C 



Due to the presence of residual corrections to finite size scaling, one has actually to extrap- 
olate the results taking the limit (Info) -1 — > in the upper part of the Fig. |8|. We observe a 
strong correction for the small sizes. However for the biggest sizes the fit seems good enough 
and we can extrapolate T C K as 

T C K = 0.56465(8) , (19) 
The estimate for the universal quantity U£ at the critical temperature is 
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U? = 0.6269(7). (20) 



This value is far away of the two dimensional Ising value U* sm9 ~ 0.611, |3tJ which is a 
strong indication that the universality class associated to the chirality order parameter is 
not of Ising type. We will verify this prediction studying now the critical exponents. 

At T = T C K the critical exponents can be determined by log-log fits. We obtain v K from 
Vy and V 2 K (Fig. |[), 7 K /z/ K from x K an d X2 (Fig- an d (3*/v K from k (see Fig. |ll|). We 
observe in all these figures a strong correction to a direct power law. It is worth noticing 
however that X% shows smaller corrections. Using only the three (four for X£ ) largest terms 
we obtain: 

v * = 0.795(20) (21) 
yyz/ K = 1.750(10) (22) 
P"/v* = 0.127(10) . (23) 

The uncertainty of T C K is included in the estimation of the errors. The large values in 
our errors come from the use of only few sizes for our fits. If we had used more, the 
exponents would change and for example v K would grow until 0.91 using all the sizes. The 
non observation of the corrections in previous studies could explain the very dispersed results 
obtained in various studies of different frustrated 2D XY models (between 0.76 to 0.90). We 
note that we have used much more statistics (due to one part to a better algorithm, and in 
other part to longer simulations) than previous works (between one or two order of magnitude 
more) which enables us to observe the finite size corrections. Moreover, we expect that the 
critical exponents written above could vary with the ratio Jij J\- It makes therefore difficult 
quantitative comparisons with other studies. Nevertheless, we can safely state that an Ising 
universality class is excluded. If we try to introduce a correction to calculate the exponents, 
for example like V* = (1 + L~^)L l l uK , we obtain uj = 1.0(2) and values for critical exponents 
fully compatible with ( |2~T| - [Z3"D . We have noticed that the exponents have a tendency to move 



away from the ferromagnetic Ising values when the size grows and thus seems to exclude a 
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crossover to the ferromagnetic Ising fixed point for larger sizes (unless it occurs at very large 
and not yet accessible size). 

The values given in ( pl| - |2"3| ) use the properties of the free energy at the critical tempera- 
ture. But an error on T C K leads to an error on the exponents, it is therefore interesting to find 
them without the help of T C K . This can be done using the whole finite size scaling region and 
the method given in f3l| . It consists to plot, for example, the susceptibility X K L^ lK / yK as 
function of U K , choosing the exponents as the curves collapse. This fit is stronger than the 
fit at the critical temperature in so far as it does not depend only on results at T C K but on a 



large region of temperature. However the errors are a little bit larger. We show in Fig. |12|-|14 
the results for three choices of r y K /v K . Obviously the result / y K /u K = 1.76 is the best one. 
With this method we arrive at r y K /i> K = 1.76(2) which is compatible with the result (|22|) 
and thus constitutes an indirect check of the critical temperature. We have verified, using 
cumulants V\ and V2 and < k >, that the results for v K and (3 K jv K are compatible with fl2~T|) 
and (|23|). 

From the scaling relation 

7 K /z, K = 2 - rj K (24) 



we obtain if = 0.25(1). The results are summarized in table |TJ. 

In conclusion the chirality order parameter does not seem to belong to the standard 
two dimensional Ising universality class. Such conclusion has already been reached in many 
studies of frustrated XY models. Nevertheless, due to the fact that the exponents could be 
J 2 / 1 J\ dependent, we cannot safely compare the results we get for the J\ — J 2 XY model 
with other frustrated XY models. However we observe that the exponents vary strongly if 
corrections are not taken into account and we suspect that this is also the case in the other 
studied models. 

If the transition belongs to a new universality class, the use of Binder's cumulant at the 
critical temperature U£ could be a new approach to track it. It should be very interesting 
to test this idea in other systems like the Villain or the triangular models. 
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E. U(l) symmetry 



We now focus on the phase transition associated to the U(l) symmetry, i.e. to the XY 
spins. The usual ferromagnetic XY model undergoes a Kosterlitz-Thouless phase transition 
driven by the unbounding of vortex-antivortex pairs. The best method to obtain reliable 
results is to use the jump of the helicity parameter defined by the answer of the system to 
a twist in one direction. The knowledge of the jump at the critical temperature allows to 
obtain T c with a good precision ||32|| . However, in our problem the presence of the chirality 
order parameter coupled with topological defects leads to a non-universal jump. This fact 
explains why this transition is usually not explored in Monte- Carlo simulations of frustrated 
XY models or, when it is, why results are not very accurate. In the following study we 
will use a method introduced in |33| using Binder's cumulant to study this transition. It 



was proved in this article that, contrary to the common belief, the Binder cumulant for 
ferromagnetic XY systems crosses for different sizes, allowing thus a rough estimate of the 
critical temperature and especially of the exponent r\ without the precise knowledge of the 
critical temperature. 

We proceed in a way similar as for the Ising order parameter, replacing k by M. We 
record the variation of U M fl26|) with the temperature for various system sizes in Fig [F[ We 
want to underline the differences between the result of U K (Fig |7|) and U M which are plotted 
with the same scale. U shows a crossing on a smaller region than U K , and at least one 



order of magnitude less than the standard XY model (see the figure 1 of |33|) We then locate 
the intersection of these curves and plot the results in the lower part of the Fig. |8|. 

Let us first consider a power law behavior at T > Tc for this system. In this case we 
have to consider a linear fit for (Info)" 1 — > 0. We observe corrections for the smallest size 
L = 20 but the others seem to converge to the temperature 

T C M = 0.56271(5) . (25) 

Secondly we consider the behavior to be exponential as in the standard XY model. In 
this case figure 2 of j33| shows that a linear fit could be wrong and that a "crossover" to 
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a different critical temperature could be observed for bigger b, i.e. greater sizes. However 
contrary to the ferromagnetic XY model, the region of crossing is so small and the different 
linear fits tend only to one critical temperature. Therefore we think that the linear fit 
works well enough. Moreover in the following we will show strong arguments in favor of the 
temperature fl2"5p. 

With the help of the critical temperature we have found an estimate of U M at T C M fitting 
the value with a law U M = + aL~ e . We obtain 

= 0.638(5) . (26) 

By log-log fit we calculate some exponents. The exponent rj could be obtained by a fit with 
X^ 1 shown in Fig ITUl We obtain here 



2 - 7] M = 1.657(5) (27) 
V M = 0.345(5) (28) 

The fit has been done discarding the two smallest sizes (L = 20 and L = 40) which show 
small corrections. This value is different from the standard XY case where rj = 0.25. No- 
tice also that it is in contradiction with the result of Monte Carlo simulations in the high 
temperature region obtained by Jose et al. || for the FFXY (which is believed to be in the 
same universality class as our model) where rj ~ 0.20 was found. To our knowledge, it seems 
one of the first times this exponent is calculated using finite size scaling. From a theoretical 
point of view the KT transition has an exponential behavior, i.e. a correlation length of the 
form £ ~ exp[S (T — T C ) _I/ ], however a power law behavior like (£ ~ (T — T c )~ u ) can not be 
excluded numerically. In the latter case the critical exponent v can be calculated with the 
cumulant V 2 M ([15]). We have obtained v M = 0.92(3). With the finite size scaling method 
we were not able to compute the exponent v in the case of an exponential behavior. 

As for the Ising order parameter, the calculation of the exponents have been done at 
the critical temperature but an error on T C A/ leads to errors on the exponents, it is thus 
interesting to find them without the help of T C M . This can be done using the same method 
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as described before. We have shown in [33] that this method is accurate enough in order to 



obtain rj whatever the type of the behavior is (power law or exponential). In Fig. |T5|-|i7| we 
show our results for three values of rj M . Obviously the value 7] M = 0.33 is the best and we 
are able to obtain: 



Tj = 0.33(2) (29) 

which is compatible with (P8[ ). Moreover this result is a strong indication that our choice 
of the critical temperature is correct. Indeed another choice leads to other non-compatible 
exponents. For example, had we chosen T C M = T C K = 0.56465 we would obtain rj RI (T = T c ) = 
0.47 which is incompatible with (f29"|). 

To sum up, we have computed for the first time the critical exponent r] M = 0.345(5) for 
the Kosterlitz-Thouless transition using the finite size scaling region in Monte Carlo simula- 
tions. We have given strong arguments that, in our range of sizes, the critical temperature 
for this transition is less than the critical temperature corresponding to the Ising-like tran- 
sition (T C M < T C K ). Note that this is in agreement with the phase diagram of the Ising-XY 
model (Fig 0), where if the two transitions never merge, we have Tkt < Tj (see also ref. 



[15]). We cannot exclude, contrary to the ferromagnetic XY model, a power law behavior 
at T > T C M which should be characterized besides the exponent r] M , by v M ~ 0.92 and 



IV. CONCLUSION 

In this paper we have investigated the critical behavior of the 2D, J\ — J 2 XY model on 
the square lattice. We have first theoretically argued that this model should be in the same 
universality class as the Ising-XY model for 2J 2 > J\. We have then carried on extensive 
Monte Carlo simulations for the particular ratio J2/J1 = 0.7. Our main conclusion is that 
this system undergoes two distinct and separate transitions. The first one is of Kosterlitz- 
Thouless type with an exponent rj = 0.345(5) different to the ferromagnetic case, whereas 
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the second one (associated to the chirality order parameter) seems to be in a non-Ising 
universality class. The temperatures of transitions and the values of the critical exponents 
are summarized in table [TTT| . It is worth mentioning that the estimate of the exponents can 
be obtained without the help of a precise determination of the critical temperatures. The 
values of the critical exponents associated to the Ising symmetry are consistent with those 
obtained in various recent works for different frustrated XY models P,|6|,|8|j9|,|T5|] . Nevertheless, 
our analysis has been done at J2/J1 = 0.7. We expect that the exponents we obtained 
could vary with the ratio J2/J1 which makes accurate comparisons difficult. Consequently, 
numerically speaking, we cannot safely state that the J\ — J2 XY model is in the same 
universality class as other models quoted above. 

The fact that two transitions exist at two different temperatures and that the critical 
exponents of the chiral order parameter transition is not of Ising type seems puzzling. How 
could we reconciliate them ? One first idea, is to use an argument by Olsson |~3] which 
explains these strong deviations from the Ising universality class by a large screening length 
(due to vortices) which prevents observing the expected Ising behavior. In our case, we 
would therefore expect for large sizes, a crossover to such behavior, i.e. for example, v grows 
to reach the value 1 for an infinite size. However no sign of this crossover has appeared for 
our largest sizes [L = 150). Obviously, such a crossover could not be excluded for much 
larger size (L ~ 1000) but seems not very plausible. Two more plausible interpretations can 
be a priori formulated. One, due to Granato et al. consists in invoking a new universality 
class for the chiral order parameter. The idea of the 3— state Potts model universality class 
has for example been recently advocated in ||. Another one is due to Knops et al. [|J. 
These authors have introduced a possible instable fixed point on the critical line PT of the 
phase diagram of the Ising- XY model (see Fig. which is able to induce strong cross over 
phenomena in the infrared limit. This conjecture has the advantage to explain the whole 
set of dispersed results found in the literature. Notably, it would explain the continuous 
variation of exponents found by Granato et al. in the Ising-XY model [^| but also the 
the J2/J1 dependence of critical exponents in the J\ — J2 XY model under consideration 
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here. To answer these questions, very high precision Monte-Carlo simulations for large size 
systems could bring some answers to this problem. Moreover new analytical developments 
are deeply needed. 
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TABLES 



NOR 



t k (1 + NOR/ 6) 



256(9) 



256(9) 



18.9(4) 



25.1(5) 



14.5(2) 



21.7(3) 



12.3(2) 



20.5(3) 



11.2(1) 



20.5(2) 



10.2(2) 



20.3(3) 



9.5(2) 



20.5(3) 



10 



8.4(1) 



22.4(2) 



15 



7.50(4) 



26.2(1) 



TABLE I. Autocorrelation times for the chirality for L = 30 in function of the number of 
over-relaxation steps NOR. 
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25.58(33) 
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1.10 6 


16.10 6 


31.66(50) 


510 5 


120 


2.10 6 


20.10 6 


40.6(10) 


510 5 


150 


3.10 6 


32. 10 6 


50.9(15) 


610 5 



TABLE II. Number of Monte Carlo steps to thermalize T up and to average T av as function 
of the size of the lattice L. t k are calculated with shorter MC runs. The last column gives the 
number of "independent" measures which are, at least, one or two orders greater than previous 
studies. 
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symmetry 
Z 2 



0.56465(8) 



0.6269(7) 



0.795(20) 



1.750(10) 



0.250(10) 1 



P/v 
0.127(10) 



U(l) 



0.56271(5) 



0.638(5) 



0.92(3) 2 



0.345(5) 



TABLE III. Summary of our results for the Ising symmetry (Z 2 ) and the XY symmetry U(l). 
1 calculated using 2 — 77 = 7/1/. 2 for a power law behavior. 
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APPENDIX A: CALCULATION OF ERRORS 



We describe here our procedure to calculate statistical errors for the different quantities. 
The first stage is to find the number of independent steps in our Monte Carlo. Indeed the 
Monte Carlo is a Markov process and therefore two consecutive steps are correlated. 

1. Autocorrelation time 

We define the autocorrelation function: 

pa) - < mm >-< Ait) 

p[t) - <A(ty>-<A(t)>* ' (A1) 

where A(t) is a thermodynamic quantity (n, \ M > X K •••)• An example is shown in figure |5| 
for the chirality k at the temperature T = T C K and for a lattice size L = 20. Nmc — 50 
millions steps of Metropolis algorithm are used after discarding 1 millions steps to thermalize 
the system. We calculate the autocorrelation time following the procedure of by r = 



Y^t=o pit) where tf is calculated in a self consistent way as pit) < 0.01 which corresponds 
to tf ~ 5t. In this case the value we get is systematically underestimated for less than 1 
percent. It is important to stop the summation because the variance of r is of the order of 
the number of summation (tf) and thus the error proportional to y/tj. Madras and Sokal 
have proposed an estimation of the error: 

At ~ y/2(2t f + l)ryN MC . ( A2) 

However this formula seems to underestimate the result. Indeed we calculate t(At) function 
of the number of MC steps. We obtain 104(5), 112(2) and 116.2(8) for one, ten and fifty 
millions respectively. Obviously the errors are underestimated (104 + 5 < 116.2 — 0.8). 
Therefore, in order to compute the errors for r we make several simulations for different 
initial configurations and use the results as independent quantities to calculate the average 
and estimate the error. 
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2. Statistical Errors 



The second step after having computed the autocorrelation time r, is to calculate the 
error on each quantities. As they depend only on a single average like the chirality or the 



susceptibility x% (HH)? the result is straightforward [35 



(A < n >f = (< > - < K > 2 ) yl ^ 2T ' Ts \ (A3) 

(A Xl f = >>)G+*^f> (A4) 

NmcHs k B ± 

where Nmc is the number of Monte Carlo steps to average, t s the number of steps between 
two measurements and N the number of lattice sites. Choosing t s = 1 we obtain (1 + 
2t)/Nmc ~ 2t/Nmc for large r while choosing r s = r we get 3t/Nmc which gives larger 
errors. 

Problems arise when quantities are a combination of different averages, for example the 
chirality \ K (FD- We could try to treat < k > and < k 2 > as independent quantities and 
estimate the error by the sum of the errors of the two quantities. However the result will 
be overestimated due to the correlation between the two elements of the sum. To solve this 
problem we can use, for example, the Jackknife procedure. We do not review this method 
here but just present the essential points we need (for a more complete review see p6|). 

In order to avoid being too abstract, we show how this method works for the susceptibility 
of the chirality (fTT|). For clarity we choose t s = 1. The application for different t s is then 
straightforward. We have to define 

_ ( = N MC < K >- Kt 

- 2( = N MC < K * > - A 
Nmc — 1 

Xt = (A7) 
k B T 

where t designs the MC step and Nmc the total number of MC steps. Our estimate for the 
susceptibility and the error will be given by: 
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AT 1 N MC 

Ax 2 ~^f^Efe-x) 2 *(l + 2r) (A9) 

iV A /C t= i 

If we save the chirality at each MC step the formulas are not difficult to apply. However we 
need a large hard disc to store the data. For example if we wanted to save the 32 millions 
steps for the simulation of the size L = 150 we would then need 72 bytes for each step 
to save the energy, the magnetization and the chirality, which implies more than two giga 
bytes. To avoid this problem, we could only save the data every r s = r but then the size of 
the file would still be more than ten millions of bytes. Moreover we would lose informations 
and therefore errors would be greater. We propose now a way which allows to obtain a good 
estimate without the problem of large storage and with minor changes in the program. 

We use a development for large n = Nmc/Q + 2r) (which is always the case in MC), 
choosing r s = 1. In this case the formula ( |A8[ - |A9| ) becomes, keeping only the dominant 
term: 

X ~( <K 2 >- <K> 2 )JL (A10) 

k B T 

A X 2 ~ (1 + 2r) [ < K 4 > - < K 2 > 2 +4 < K > 2 (< K 2 > - < K > 2 ) 

N MC 

-4 < K > (< K 3 > - < K 2 >< K >)](y^-) 2 (All) 

KbT 

The chirality conserves its initial form while the error is the sum of the two errors (of < k 2 > 
and < k > 2 ) subtracted by the third term which representing the correlation between them. 
We note that this procedure induces a small change in the program: we have only to save 
in the histogram < k 3 > plus the values of < k > and < k 2 >. To test our formula we have 
computed the susceptibility associated to the chirality ([□]) and its errors calculated by three 
ways. We perform the simulation with a lattice size L = 20 with 4 steps of over relaxation 
algorithm between each Monte Carlo, for one million steps. In this case the autocorrelation 
time is about 8 (see table [Tip. The first method consists in saving at each step the energy, 
the magnetization and the chirality, the second in saving the data only at each r steps, while 
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the third in using the approximate formula (|A11|) . We obtain, respectively: x K = 11-22(8), 
11.23(10), 11.22(8). The three methods give compatible results but the third one gives the 
best estimate with the smallest size of storage (some hundred thousands bytes). 

We give hereafter the results of our calculation for the binder parameters ( |T6| - |T7| ) and 
the cumulant V x and V 2 (f§0): 

A f/ 2 ~ ( \+ 2r) [4 < K A > 3 -4 < K 2 >< K A >< K 6 >+8< K 2 > 2 < K S > 

Nmc 

-<k 2 > 2 < k 4 > 2 } (A12) 
AV 2 „ (l + 2r) <kE > 2 <k 2 E 2 > _ 2 < k 2 E > + <k 2 > 
1 N M c < « > 2 < kE > 2 < kE >< k> < k> 2 
+ < E 2 > - < E > 2 

-2^^( < ^ 2> - <^>)] (A13) 

<K> ^ <KE> < K> Ji V; 

AV2 = AVi substituing n by k 2 (A14) 



3. Systematic Errors 

In addition to statistical errors, we have to take care of systematic errors. There are 
essentially of two kinds: one due to the correlation between the random number and one 
due to the use of the histogram technic. 

The first one appears when we use a bad random generator. In this case the period of the 
random numbers could be very small and could thus introduce correlations between data. 
One example is the linear congruential generator used by many physicist for Monte Carlo 
simulations! For certain choices of the initial parameter, the period could be very small (less 
than 2000) and therefore could induce systematic errors. We use in this work a random 
generator with a period of more than 100 millions found in Numerical Recipes (rani) ||36|| . 

A second source of systematic error comes from the histogram technic and the difference 
between the temperature of simulation and the temperature where the quantities are com- 
puted |3^J3^|. In our simulation this difference is kept less than 0.005 in order to minimize 



this effect. 
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FIGURE CAPTIONS 



Fig. [I]: Ground state of the J\ — J 2 XY model for 2 J 2 > J\ 

Fig. 0: Phase diagram of the Ising-XY model. Solid and dotted lines indicate continuous 
and first-order transitions respectively. 

Fig. [| Real autocorrelation time for the standard Metropolis algorithm (circle) and in 
combination with the over- relaxation algorithm (square). 

Fig. IJ Phase diagram for the J\ — J 2 model. For 2 J 2 < J\ we find the normal Kosterlitz- 
Thouless transition. Lines are just guides for the eyes. Our study is done for J%/ J\ = 0.7 
(black circle). 

Fig. |5|: Autocorrelation p(t) for the chirality k at T = T C K . The lattice size is L = 20 and 
the number of MC is Nmc — 50 millions. The estimated r = 116(2) is shown by an arrow. 

Fig. |6]: Binder's parameter U K for the Ising order parameter function of the temperature 
for various sizes L. The arrow shows the temperature of simulation T s = 0.565. 

Fig. ^ Binder's parameter U M for the U(l) order parameter function of the temperature 
for various sizes L. The arrow shows the temperature of simulation T s = 0.565. the scales 
is similar to those of Fig. |6| 

Fig. g: Crossing T plotted vs inverse logarithm of the scale factor b = L'/L. The upper 
part of the figure corresponds to U K while the lower part to U M . In the last case the size 
L = 60 is not shown for clarity. We obtain T C K = 0.56465(8) and T C M = 0.56271(5) with a 
linear fit (see text for comments). 

Fig. |9|: Values of V" and V£ function of L in a log-log scale at T C K . The value of the 
slopes gives \jv K . We observe strong corrections for small sizes. Only the three largest sizes 
are used for the fits. When not shown, the estimated statistical errors are smaller than the 
symbol. 

Fig. [H]: Values of \ K an d X2 function of L in a log-log scale at T C K and x¥ a ^ T^ 1 . The 
value of the slopes gives 7/1/ = 2 — rj. We observe strong corrections for the small sizes 
for x K - Only the three largest sizes are used for the fit for x K while only the smallest sizes 
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L = 20 and L = 40 are discarded for the fits for x% an d x¥ ■ When not shown, the estimated 
statistical errors are smaller than the symbol. 

Fig. [Tl]: Values of < k > as function of L in a log-log scale at T C K . The value of the 
slopes gives (3 K jv K . We observe strong corrections for the small sizes. Only the three largest 
sizes are used for the fits. When not shown, the estimated statistical errors are smaller than 
the symbol. 

Fig. U: x K L~ r/uK function of U K with 7 K /z/ K = 1.79 for the sizes L = 60, 80, 100, 120 
and 150. The curves do not collapse in one curve. 

Fig. p): x K L^ K/uK as function of U K with 7 K /z/ K = 1.76 for the sizes L = 60, 80, 100, 120 
and 150. The curves collapse in one curve. 

Fig. 0: x K L^ K/uK function of U K with ^jv K = 1.73 for the sizes L = 60, 80, 100, 120 
and 150. The curves do not collapse in one curve. 

Fig. ||: x M L 2 ~ vM function of U M with rj M = 0.31 for the sizes L = 60, 80, 100, 120 
and 150. The curves do not collapse in one curve. 

Fig. [TB|: X M L 2 -^' function of U M with r) M = 0.33 for the sizes L = 60, 80, 100, 120 
and 150. The curves collapse in one curve. 

Fig. [T7|: x M L 2 ^ M function of U M with rj M = 0.35 for the sizes L = 60, 80, 100, 120 
and 150. The curves do not collapse in one curve. 
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